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Abstract.  The  objective  of  this  study  is  to  broaden  the  fundamental  understanding  of  the  emerging  field  of  microfluidics 
especially  in  a  long  channel.  The  quasi  gasdynamic  (QGD)  equations,  originally  developed  on  the  basis  of  a  kinetical  model 
are  used  for  numerical  and  analytical  simulation.  A  two-dimensional  analysis  of  the  QGD  equations  with  a  first  order  slip 
velocity  boundary  conditions  demonstrates  that  both  compressibility  and  rarefied  effects  are  present  in  long  microchannels. 
Analytical  solutions  for  the  pressure  and  the  velocity  profiles  are  derived  from  the  quasi  gasdynamic  equations  by  undertaking 
perturbation  expansions  according  to  a  small  parameter  £  (the  height-to-length  ratio  of  the  channel)  and  using  the  isothermal 
flow  assumption.  The  deduced  expression  for  the  mass  flow  rate  is  similar  to  the  analytical  expression  obtained  from  the 
Navier- Stokes  equations  with  a  second  order  slip  boundary  condition  and  gives  results  in  agreement  with  the  measurements. 
The  effects  of  the  rarefaction  and  of  the  compressibility  on  pressure  distributions  are  analyzed.  The  analytical  expression  of 
the  pressure  predicts  accurately  the  measured  pressure  distribution.  The  Knudsen  numbers  calculated  at  the  exit  of  the  channel 
and  based  on  the  channel  height  vary  from  10-3  to  0.4.  The  comparisons  of  analytical  and  numerical  solutions  confirm  the 
validity  of  the  analytical  approach. 


INTRODUCTION 

A  systematic  research  effort  in  micro  mechanics  devices  started  in  the  late  1980’s.  The  characteristic  length  scale 
of  these  devices  is  less  than  1  mm.  Then,  even  at  the  atmospheric  conditions,  the  ratio  of  the  mean-free  path  to 
the  characteristic  dimension  can  not  be  neglected  and  in  the  flow  dynamic  associated  with  MEMS  the  rarefied  gas 
phenomena  become  apparent.  In  the  present  study  an  analytical  and  numerical  investigation  of  slightly  rarefied  gaseous 
flows  through  long  microchannels  is  undertaken.  Analytical  expressions  of  the  streamwise  mass  flow  rate,  pressure 
and  velocity  profiles  are  derived  from  the  quasi  gasdynamic  (QGD)  equations  [1]  through  perturbation  expansions 
according  to  £  (the  height-to-length  ratio  of  the  channel)  under  the  isothermal  flow  assumption. 

We  obtain  a  mass  flow  rate  expression  similar  to  that  previously  derived  from  the  Navier-Stokes  (NS)  equations  as¬ 
sociated  to  a  second  order  slip  boundary  conditions.  This  expression  gives  results  in  agreement  with  the  measurements. 
Regarding  this  point,  we  must  especially  mentioned  a  study  [2]  based  on  the  quasi  hydrodynamic  (QHD)  equations 
written  out  for  a  plane  stationary  isothermal  gas  flow  in  a  long  channel;  since  the  QHD  and  QGD  systems  are  both 
founded  on  the  same  basic  concept,  even  if  involving  different  closures  of  the  conservation  equation  system  [1],  thus 
an  analytical  mass  flow  rate  formulation  was  given  in  [2]  very  close  to  that  presented  in  this  paper.  Furthermore  we 
derive  from  our  perturbation  method  the  pressure  and  the  velocity  profiles  at  £  zero  order.  The  analytical  expression 
of  the  pressure  predicts  accurately  the  measured  pressure  distribution  [3].  The  effect  of  the  rarefaction  and  of  the  com¬ 
pressibility  on  the  pressure  and  on  the  velocity  distributions  are  pointed  out  and  discussed.  Otherwise  the  numerical 
calculations  based  on  the  full  system  of  the  QGD  equations  were  carried  out  for  different  sizes  of  the  microchannels 
and  for  various  gases.  The  agreement  between  analytical  and  numerical  results  confirms  the  validity  of  the  analytical 
approach. 


QGD  EQUATIONS:  SIMPLIFIED  ISOTHERMAL  FORM 

The  QGD  equations  are  a  gasdynamic  system  closed  with  a  special  choice  of  mass  flux  vector,  shear-stress  tensor 
and  heat  flux  vector,  determined  by  means  of  time- space  averaging  (instead  of  space  averaging  for  the  Navier-Stokes 
equations)  for  the  gasdynamic  quantities  p,  Ui  and  p  [  1].  The  general  form  of  the  QGD  system  is  presented  in  [1], 
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[4].  The  2-D  flow  is  considered  to  be  isothermal,  therefore  the  energy  equation  is  not  taken  into  account  in  the  present 
analysis. 

The  variables  are  non-dimensioned  in  the  following  manner:  the  streamwise  coordinate  v  by  the  channel  length  L; 
the  wall-normal  coordinate  y,  by  the  channel  height  //,  the  velocities  u  and  v  are  normalized  by  the  velocity  u0  on 
the  axis  at  the  channel  exit  (the  subscript  "o"  refers  to  the  outlet  conditions);  the  density  and  the  viscosity  coefficient 
by  the  appropriate  outlet  values  p0  and  p0.  The  channel  height-to-length  ratio  £  will  be  small  compared  to  unity: 
£  =  H /L,£  <C  1.  The  choice  of  the  pressure  normalization  is  based  on  the  following  reason:  the  pressure  changes 
due  to  inertia  indeed  are  very  small  at  small  Mach  numbers  (that  is  the  case  in  microflows  [5]).  Thus  in  a  long 
microchannel  these  pressure  changes  along  the  channel  are  mostly  due  to  viscous  effects.  Therefore  we  suppose  that 
they  are  proportional  to  the  viscous  forces:  ^  ~  >  that  leads  to  choose  as  the  referenced  pressure  u0Lp0/H 2  for 

the  non-dimensionalisation.  Moreover,  in  the  various  experiments  investigated  hereafter  the  pressure  change  is  of  the 
same  order  as  the  pressure  itself:  so,  the  non-dimensional  pressure  is  guaranteed  to  be  of  zero  order  in  £  as  the  others 
non-dimensional  parameters  are.  Now  we  expand  the  velocities,  the  density  and  the  pressure  in  powers  of  £  up  to  the 
first  order: 


u  =  uo  +  £u\ ,  v  =  £vi,  p=Po  +  epi,  P  =  Po  +  £P1j 


(1) 


and  we  substitute  the  expressions  (1)  into  the  non-dimensional  system  as  it  was  fulfilled  in  similar  analysis  of  the 
NS  equations  [6].  Disregarding  the  flows  of  0(  1/e)  Reynolds  numbers  (where  Re  =  p0u0H/ju0,  evaluated  in  the 
outlet  conditions)  analyzed  in  the  boundary  layer  theory  [7],  the  remaining  flow  regimes  characterized  by  Re  ~  0(1) 
and  0(e)  are  regimes  classified  as  microflows  and  may  be  appropriately  treated  by  the  current  perturbation  analysis. 
Multiplying  the  y-momentum  equation  by  the  factor  £2Re  and  keeping  the  £  zero  order  terms  we  obtain  (omitting  from 
now  the  zero  subscript)  dp/dy  =  0,  thus  the  pressure  is  expected  to  be  uniform  across  any  section  of  the  channel.  We 
use  this  result  and  the  isothermal  flow  assumption  (so  the  temperature  and  viscosity  coefficient  are  constant),  and  we 
neglect  the  terms  of  £/Re,  1  and  1  /Re  order  against  the  term  of  1  /(sRe)  order.  Then  multiplying  the  v-momentum 
equation  by  £Re  we  obtain  at  zero  order 


dp  „  d2u 
dx  ^  dy2 


(2) 


We  can  remark  that  the  momentum  equations  are  the  same  as  those  obtained  with  the  NS  system  [6],  [8]. 

Then,  the  boundary  conditions  for  the  equation  (2)  are  the  following:  the  symmetry  condition  on  the  axis  and  the 
slip  velocity  boundary  condition  on  the  solid  wall,  that  in  the  isothermal  case  reads  [9] 


(3) 


where  vm  =  \flMTw  is  the  most  probable  molecular  velocity  at  the  surface  temperature  Tw\  op  is  the  slip  velocity 
coefficient.  The  value  of  the  slip  velocity  coefficient  op  suggested  by  Maxwell  is  o^  =  0.5v^(2  —  a) /a,  where  a 
is  the  part  of  the  molecules  reflected  diffusively.  In  the  case  of  the  full  accommodation  the  coefficient  is  equal  to 
0.886.  But  in  this  article,  unless  advised  otherwise,  we  will  use  for  the  QGD  equations  the  velocity  slip  coefficient 
off  =  1.012,  given  by  Kogan  [10]  under  the  full  accommodation  assumption.  This  value  is  close  to  the  one  proposed 
by  Sharipov  [9]  o^  =  1.0  in  the  case  of  the  diffuse  gas-surface  interaction. 

Furthermore  the  mean  free  path,  as  a  function  of  macroscopic  parameters,  is  usually  written  as  A  =  k^p/ pVWf.  In 
this  equation  the  coefficient  k %  depends  on  the  molecular  interaction  model.  For  the  HS  model  the  Chapman  formula 
[11]  leads  to  k ^  =  \Jn/2.  For  the  present  QGD  equations  we  will  use,  unless  advised  otherwise,  =  A(co)  = 

2(7~15v^2C°)  ’  exPressi°n  deduced  by  Bird  [12]  for  the  VHS  model;  the  coefficient  A(co)  depends  only  on  the 
type  of  the  gases.  The  slip  velocity  is  related  to  the  mean  free  path  and  to  the  velocity  gradient  through  the  expression: 


y/2  f  du\  f  du\  f  du\  y/2 

Us  =  op—Ay-^yj  =  KslipA  ^  —  J  ,  or  in  non-dimensional  form  us  =  KsUpKn  j  ,  KsUp  =  op  —  ,  (4) 

where  Kn  =  A  /H  is  the  local  Knudsen  number  and  Ksup  is  a  non-dimensional  coefficient. 

Using  the  symmetry  condition  and  the  slip-flow  boundary  condition  (4)  the  simplified  v-momentum  equation  (2) 
can  be  integrated  twice  with  respect  to  y,  leading  to: 

U (x. ,59  =  -  2.  ( ] I  -  4y2  +  AKslipKn)  .  (5 ) 
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The  same  form  of  the  velocity  distribution  was  obtained  using  a  first  order  boundary  condition  associated  to  the  Navier- 
Stokes  equations  in  [6],  [8].  Note  that  the  streamwise  velocity  depends  explicitly  upon  the  wall-normal  (y)  direction 
and  also  on  the  streamwise  (x)  direction  via  the  local  Knudsen  number. 

We  can  now  calculate  the  mass  flow  through  the  channel  for  given  inlet  and  outlet  pressures.  Neglecting  the  terms 
of  order  equal  or  superior  to  0(e)  in  the  non-dimensional  form  of  the  QGD  mass  flux,  we  obtain: 


Jx  =  pu  - 


£  jl  dp 
Re  p  dx 


(6) 


In  comparison  with  the  mass  flux  given  by  the  NS  equations  the  QGD  mass  flux  contains  the  additional  term 
proportional  to  the  streamwise  pressure  gradient.  This  term  is  kept  at  zero  order  assuming  Re  ~  0(e).  Integrating 
the  mass  flux  (6)  across  the  channel,  using  the  streamwise  velocity  (5)  and  the  state  equations,  we  obtain: 


J=- 


1  dp  yMa2 
24 pT  dx  eRe 


2p+\2Ksupkx 


■sfyMa 


Vt  +  24 


fi2f 


yMa2  p 


(7) 


Integrating  the  mass  flow  rate  (7)  in  the  streamwise  direction  from  0  to  1  and  expressing  the  local  Knudsen  number  as 
Kn  =  Kn0p0/ p ,  then  using  the  mass  flow  conservation  along  the  channel,  we  obtain  finally  a  mass  flow  rate  expression 
depending  on  the  inlet  and  outlet  pressure  &  =  pi/ pQ\ 


pi  yMa2 
24  pT  eRe 


24  > 

*2  -  1  +  12 KslipKn0  (£»-!)  +  -^Knl ln&> 


(8) 


This  relationship  presents  in  fact  the  non-dimensional  flow  rate  per  width  unit.  Then  the  effects  of  the  slip  condition 
and  of  the  QGD  additional  term  in  the  mass  flux  appear  more  clearly  dividing  (8)  by  the  mass  flow  rate  deduced  from 
the  NS  equations  with  no-slip  condition  (J™fip): 


Jqgd/JZUP  =  1  +  12  KsUpKn0 


24  , 

+  T2Kno 


In  3^ 

^2-l’ 


jnoslip 

JNS 


24/17’  eRe  v  ’ 


(9) 


Moreover  using  the  first  order  slip  boundary  condition  with  the  NS  equations,  the  mass  flow  expressed  by  [6],  [8],  [13] 
takes  the  form 


fW/jnoslip  =  !  +  l2KslipKno_L_'  (10) 

The  last  term  of  expression  (9)  (in  Kn 2  order)  is  the  additional  term  introduced  by  the  QGD  model  in  respect  to  the  NS 
model  (10),  when  both  the  models  are  used  with  the  same  first  order  boundary  condition.  As  said  in  the  introduction, 
this  term  is  similar  to  that  derived  from  QHD  model  in  [2] :  the  two  expressions  differ  only  by  the  Schmidt  number 
which  appears  in  denominator  of  the  Kn2  term  of  QHD  model.  Note  that  in  [2]  comparing  their  results  with  the 
measurements  the  authors  found  an  agreement  in  the  0.1  —  0.5  Kn  range. 

Various  expressions  of  second  order  slip  boundary  conditions  were  suggested  by  numerous  authors  [8],  [13]  to 
supplement  the  NS  equations:  so  in  [8]  a  second  order  boundary  condition  leads  to  the  following  mass  flow  rate 


jslipl  /  jnoslip 
JNSb  /JNS 


1  + 12  KsnpKn0- 


9  4-1 


-12- 


2-  a 


a 


-Kn 


In; 


°  &2  —  1 


(11) 


The  authors  [8]  explained  that  the  effect  of  the  second  order  correction  is  to  reduce  the  increase  in  mass  flow  rate  due 
to  the  first  order  slip.  On  the  contrary,  implementing  in  the  NS  equations  another  second  order  boundary  conditions 
[14],  other  authors  [13]  obtained  a  correction  increasing  the  mass  flow  rate: 


jslipl  j  jnoslip 
JNSC  /JNS 


1  + 12  KsnpKn0 


1 

^+1 


+  21  Kn2 


ln^ 

^2-l’ 


(12) 


and  this  expression  seems  to  lead  to  results  in  agreement  with  the  measurements  of  the  authors. 
In  Fig.  1  a)  are  compared  for  helium  flows: 


•  Three  results  obtained  in  the  Navier-Stokes  theoretical  frame,  utilizing  respectively  the  expression  (10)  derived 
from  the  first  order  formula,  the  mass  flow  rate  (11)  coming  from  a  second  order  slip  condition  [8]  and  the 
expression  (12)  constructed  in  [13]  from  another  second  order  slip  condition  [14]. 
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•  The  results  obtained  from  the  present  analytical  QGD  treatment  coupled  to  the  first  order  slip  condition  which 
leads  to  expression  (9). 

•  The  corresponding  experimental  results  [6]. 

These  calculations  and  measurements  concern  helium  streams  through  a  1.33/im  deep  rectangular  microchannel 
and  are  obtained  for  the  atmospheric  outlet  pressure  and  Kn0  =  0.165.  In  order  to  obtain  significant  comparisons  we 
have  firstly  calculated  the  various  theoretical  mass  flow  rates  with  the  same  "slip  parameters":  i.e.  kx  =  y^/2  and 

op  =  ^  ,  and  in  any  case  a  complete  diffuse  reflection.  In  these  conditions  the  QGD  approach  agrees  with  the 

measurements  slightly  better  than  the  other  approaches  (Fig.  1  a).  Then  we  have  plotted  the  QGD  results  calculated 
with  our  usual  "slip  parameters"  (kx  =A(g> ),  o^  =  1.012)  admitting  again  a  complete  diffuse  reflection.  One  can  see 
in  Fig.  1  a,  that  the  improvement  is  at  the  least  of  the  same  order  when  changing  the  slip  parameters  (kx ,  op)  as  when 
changing  the  Knudsen  number  order  of  the  continuum  approach.  Furthermore,  the  prevalent  governing  effects  of  the 
solid  surface  parameters  appears  clearly  if  noting  that  kx  disappears  from  all  the  coefficients  (Kn0Ksup  and  Kn0/kx ) 
in  the  previous  mass  flow  rate  expressions  (9)-(12),  while  op  remains  explicitly.  So,  the  reflection  process  at  the  wall 
exerts  here  a  direct  influence  while  the  intermolecular  forces  act  only  through  the  viscosity  coefficient,  whatever  the 
interaction  model  used  in  the  gas. 


STREAMWISE  PRESSURE  DISTRIBUTION  AND  VELOCITY  PROFILES 

By  integrating  the  mass  flow  rate  (7)  in  the  streamwise  direction  from  0  to  x  we  obtain  an  implicit  relation  for  the 
pressure  distribution  along  the  channel: 

(«j)  +  «^.p|  |)  +24^*2fl"^  =  -24^f|s-  <13) 

Accounting  for  p(x)  <  pi,  we  expand  the  pressure  logarithmic  function  in  the  Taylor  series  of  pressure  difference.  If 
we  truncate  the  above  equation  to  retain  only  the  terms  of  the  first  order  in  pressure  difference,  we  obtain  a  simple 
explicit  expression  for  the  pressure  distribution: 


^  \  (^  +  6 KsUpKn0)2 -xF(^)+ ^ ^  +  6KsUpKn0)  -  (^>KsUpKn0  1  12  Kn° 


Po 


&  k\ 


24  Kni 


F(0>)  =  -  1  +  12  KslipKno(0» +  — 

A 


(14) 


Using  the  expansion  of  the  logarithmic  function  up  to  the  second  order  terms,  leads  to  a  more  accurate  explicit 
expression  of  the  pressure.  The  analytical  expression  obtained  in  [6]  from  the  NS  equations  with  the  first  order  slip 
boundary  conditions  gives: 


P(x)/po  =  +  6K^~K^)2  ~x  (&1  -  1  +  UK^KnM?  -  1))  -  6 KslipKn0.  (15) 

Comparing  the  last  two  expressions,  it  appears  that,  although  both  parabolic,  the  pressure  distributions  (14)  and  (15) 
present  different  curvatures  along  x.  Furthermore  if  expanding  the  square  root  in  Kn  power  only  the  Kn  zero  order  terms 
would  be  the  same  in  (14)  and  (15).  As  previously  noted  for  the  mass  flow  rates,  the  pressure  profiles  are  independent 
of  kx  -  On  the  other  hand,  in  Fig.  1  b  the  analytical  pressure  distribution  (14)  is  compared  with  pressure  measurements 
[3],  corresponding  to  a  1.2 fim  deep  channel  charged  with  nitrogen  and  exhausting  to  atmospheric  conditions:  the 
predicted  pressure  distribution  is  non  linear  and  coincides  very  well  with  the  measured  values. 

It  may  easily  be  proven  from  (14)  that  the  curvature  of  the  pressure  distribution  is  negative,  since 

dp(x)  _  F{0>) 


(&  +  6KslipKn0) 


dx 


1(^  +  6  KslipKn0)2-xF(^)+(§^ 
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Mass  flow  (kg/s) 


FIGURE  1.  a)  Helium  mass  flow  in  1.33  micron  height  channel  Kn0  =  0.165.  The  circles  and  triangles  are  the  experimental 
measurements  of  [6].  The  solid  line  with  the  circles  presents  the  QGD  approach  (9),  the  solid  curve  is  the  NS  solution  with  the 
first  order  slip  condition  (10),  the  dotted  curve  is  the  NS  solution  with  the  second  order  slip  condition  proposed  by  [8]  (11),  the 
dash-dotted  curve  is  the  NS  solution  with  the  second  order  slip  boundary  condition  proposed  by  [13]  (12).  For  all  the  expressions 
mentioned  above  the  Chapman  coefficient  =  y/n/2  and  the  Maxwell  value  for  op  are  used.  The  dashed  curve  presents  the  QGD 
model  with  usual  "slip  parameters"  =  A(co ),  op  =  1.012  (9).  b)  Pressure  distribution  in  a  long  microchannel.  The  solid  lines 
represent  the  analytical  solutions  (14)  for  three  different  pressure  ratios  &  =  2.02,2.36,2.70,  the  symbols  are  the  data  from  [3]. 


The  limit  of  the  pressure  gradient  when  Kn  —>  0  tends  toward  the  expression  (1  —  S?2') / y/ 2?2  —  x(&2  —  l)/2  which 
depends  on  the  distance  x.  That  means  that  the  compressibility  effects  remain  when  the  rarefaction  effects  disappear 
and  that,  for  the  least,  a  part  of  the  nonlinearity  of  the  pressure  distribution  is  due  to  the  change  in  the  density  of  the 
gas  in  the  channel  and  not  to  the  rarefaction  effects.  Otherwise,  it  is  clear,  from  Fig.  2  a,  that  presents  the  deviation 
of  the  streamwise  pressure  from  linear  pressure  drop,  and  from  various  numerical  tests  that  for  the  large  Knudsen 
numbers  (in  the  frame  of  the  continuum  approach)  the  pressure  distribution  becomes,  at  last,  linear.  So,  it  could  be 
deduced  that  the  rarefaction  effects  serves  generally  to  decrease  the  curvature  in  the  pressure  distribution  caused  by 
the  compressibility  effect. 

Finally,  let  us  note  that  expression  (16)  allows  to  obtain  explicit  velocity  profiles  (putting  (16)  in  (5))  which, 
according  to  Fig.  2  b,  agree  to  the  numerical  calculation. 


CONCLUSION 

We  have  focused  our  interest  on  the  isothermal  continuum  approaches  in  the  slip  regime  using  the  quasi  gasdynamic 
equations  to  delimit  their  validity  domain.  From  the  QGD  system,  we  have  derived  a  mass  flow  rate  expression  and 
analytical  profiles  of  pressures  and  velocities,  using  a  small  perturbation  treatment.  Up  to  Knudsen  numbers  close  to 
0.2  the  QGD  approach  gives  results  fitting  slightly  better  the  measurements  than  the  NS  approaches  (employed  with 
first  or  second  order  slip  boundary  conditions).  Moreover  in  this  Knudsen  range  our  analytical  profiles  agree  perfectly 
with  the  "exact"  QGD  numerical  results.  In  this  range  three  main  results  may  be  pointed  out: 

•  The  QGD  mass  flow  rate  expression  analytically  deduced  from  QGD  equations  is  similar  to  those  obtained  from 
the  Navier-Stokes  equations  associated  to  a  second  order  slip  boundary  conditions. 

•  The  rarefaction  effects  reduce  the  curvature  of  the  parabolic  streamwise  pressure  profiles  usually  found  in  this 
Knudsen  number  range:  when  these  rarefaction  effects  disappear  (Kn0  — ►  0),  compressibility  effects  remain  and 
insure  a  maximal  parabolic  curvature  to  the  pressure  profiles. 

•  A  change  in  the  "slip  parameters"  retained  for  the  calculations  (especially  a  change  in  the  choice  between  the 
slip  velocity  coefficients  <7^  and  erf)  leads  to  result  variations,  for  the  least  as  important  as  those  observed  when 
changing  the  Knudsen  number  order  (first  or  second)  of  the  continuum  approaches. 
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FIGURE  2.  a)  Deviation  from  the  linear  pressure  drop  for  helium  flow  &  —  2.36  for  the  outlet  Knudsen  numbers  Kn0  = 
0. 164 -y  0.4.  b)  Velocity  variation  on  the  channel  axis  for  nitrogen  flow,  the  solid  line  corresponds  to  the  analytical  solution,  the 
dashed  line  is  the  numerical  solution. 


Calculations  about  more  numerous  experimental  results  and  further  investigations  on  the  slip  parameters  (accom¬ 
modation  coefficients,  reflection  law)  should  clarify  many  questions  still  arrising  in  the  microflow  field. 
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